\documentclass{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}
%\usepackage{ctex}
\usepackage{listings}
\usepackage{subfigure}


% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}

\usepackage{xcolor}
\usepackage{fontspec} 
\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{comment}{rgb}{0.56,0.64,0.68}

\newfontfamily\monaco{Monaco}
\lstset {
aboveskip=3mm,
belowskip=3mm,
showstringspaces=false,       % underline spaces within strings
columns=flexible,
framerule=1pt,
rulecolor=\color{gray!35},
backgroundcolor=\color{gray!5},
basicstyle={\small\monaco},           % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\tiny\monaco\color{gray},  % the style that is used for the line-numbers
numbersep=5pt,                  % how far the line-numbers are from the code
commentstyle=\color{comment},
keywordstyle=\color{blue},
stringstyle=\color{dkgreen},
tabsize=2,                      % sets default tabsize to 2 spaces
captionpos=b,                   % sets the caption-position to bottom
breaklines=true,                % sets automatic line breaking
breakatwhitespace=false,        % sets if automatic breaks should only happen at whitespace
escapeinside={\%*}{*)},            % if add LaTeX within your code
morekeywords={*,...}               % if add more keywords to the set
}


\begin{document}
\title{Homework \#4}
\pagestyle{fancy}
\lhead{Name Li HuiTeng 3180102114}
\chead{ numPDE\#5}
\rhead{Date 21.06.13}


\tableofcontents

\newpage

\section{10.8 Determine the matrix form of Crank-Nicolson.}
\begin{proof}[Proof]
	This follows directly from (10.15).
	Since $AU^n_i=\frac{\nu}{h^2}(U^n_{i-1}-2U^n_{i}+U^n_{i+1})$ , $\forall i=2,\cdots,m-1$we have
	\begin{align*}
		U^{n+1}_{i}-U^n_{i} & =\frac{k\nu}{2h^2}(U^n_{i-1}-2U^n_{i}+U^n_{i+1}+U^{n+1}_{i-1}-2U^{n+1}_{i}+U^{n+1}_{i+1}), \\
		U^{n+1}_{i}-U^n_{i} & =\frac{k}{2}(AU^n_{i}+AU^{n+1}_{i}),
	\end{align*}
	$\forall i=2,\cdots,m-1$.
	Also, considering
	\begin{align*}
		\frac{U^{n+1}_{1}-U^{n}_{1}}{k} & =AU^n_1+\frac{\nu}{h^2}g_0(t), \\
		\frac{U^{n+1}_{m}-U^{n}_{m}}{k} & =AU^n_m+\frac{\nu}{h^2}g_1(t),
	\end{align*}
	we have
	\begin{align*}
		(I-\frac{k}{2}A)U^{n+1}=(I+\frac{k}{2}A)U^{n}+\frac{k\nu}{2h^2}(g_0(t_n)+g_0(t_{n+1})\quad 0 \quad \cdots g_1(t_n)+g_1(t_{n+1}))^T,
	\end{align*}
	which completes the proof.
\end{proof}


\section{10.11 Show Crank-Nicolson is 2rd accurate in time and space.}
\begin{proof}[Proof]
	We have
	\begin{align*}
		                        & \tau(x,t)                                                                                                       \\
		                        & =\frac{u(x,t+k)-u(x,t)}{k}-\frac{\nu}{2h^2}(u(x-h,t)-2u(x,t)+u(x+h,t)                                           \\
		                        & +u(x-h,t+k)-2u(x,t+k)+u(x+h,t+k))                                                                               \\
		                        & =u_t+\frac{1}{2}ku_{tt}+\frac{1}{6}k^2u_{ttt}+O(k^4)-\frac{\nu}{2}(u_{xx}+\frac{1}{12}h^2u_{xxxx}+u_{xx}(x,t+k) \\
		                        & +\frac{1}{12}h^2u_{xxxx}(x,t+k)+O(h^4))                                                                         \\
		                        & =u_t+\frac{1}{2}ku_{tt}-\frac{\nu}{2}(u_{xx}+u_{xx}(x,t)+ku_{xxt}(x,t)+O(k^2+h^2))+O(k^2+h^2)                   \\
		(u_{t}=\nu u_{xx})\quad & =\nu u_{xx}+\frac{\nu^2}{2}ku_{xxxx}-\frac{\nu}{2}(2u_{xx}+k\nu u_{xxxx}(x,t))+O(k^2+h^2)                       \\
		                        & =O(k^2+h^2),
	\end{align*}
	which completes the proof.
\end{proof}

\section{10.25 Prove the necessity part of Lax Equivalence Theorem.}
\begin{proof}[Proof]
	Consider the homogeneous boundary condition and let the exact solution be u = 0. Then we have
	$U^{n+1}=B(k)^{n}U^0$ and $E^{n+1}=U^{n+1}=B(k)^{n}U^{0}$. Due to the convergence of MOL, for any T,
	we have
	$U^{n}$ inclines to zero-vector as $k \to 0$ with $(n-1)k=T$ and with $U^{0}$ inclining to zero-vector.
	Then take $f_{k,n}(v)=B(k)^{n-1}v$. In fact $f_{k,n}$ is continuous at v = 0 when $k(n-1)$ is bounded.
	This is because as k tends to 0, f maps any vector inclining to zero to a sufficiently small vector inclining to 0.
	Also, combining with the linearity of f, we conclude the fact that f is a bounded functional, which means
	$||B(k)^{n-1}||$ is bounded for any (n,k) pair satisfying $(n-1)k \leq T $ which leads to $nk \leq T$ and completes
	the proof.
\end{proof}

\section{10.28 Show the modulus of Crank-Nicolson's ampliﬁcation factor is never greater than 1 for any k,h.}
This follows directly from (10.26). We have
\begin{align*}
	g(\xi)(1+\frac{2\nu k}{h^2}sin^2(\frac{\xi h}{2})) & =1-\frac{2\nu k}{h^2}sin^2(\frac{\xi h}{2})                                                     \\
	g(\xi)                                             & =\frac{1-\frac{2\nu k}{h^2}sin^2(\frac{\xi h}{2})}{1+\frac{2\nu k}{h^2}sin^2(\frac{\xi h}{2})}.
\end{align*}
To guarantee $|g(\xi)| \leq 1$, we take a=$\frac{2\nu k}{h^2}sin^2(\frac{\xi h}{2})$. Since
$g(\xi)=\frac{1-a}{1+a}(a\geq0)$, $g(\xi) \leq 1 \Rightarrow 1-a\leq1+a \Rightarrow a\geq 0$ and
$g(\xi) \geq -1 \Rightarrow -1\leq1$, which completes the proof.

\section{11.26 Show Beam-Warming is 2rd accurate in time and space.}
\begin{proof}[Proof]
	For $a\geq0$, we have
	\begin{align*}
		                     & \tau(x,t)                                                                                      \\
		                     & =\frac{u(x,t+k)-u(x,t)}{k}+\frac{\mu}{2k}(-4u(x-h,t)+3u(x,t)+u(x-2h,t))                        \\
		                     & -\frac{\mu^2}{2k}(u(x,t)-2u(x-h,t)+u(x-2h,t))                                                  \\
		                     & =u_t+\frac{1}{2}ku_{tt}+O(k^2)+\frac{a}{2h}(4hu_x-4(-h)^2u_{xx}/2+O(h^3)-2hu_x+(2h)^2u_{xx}/2) \\
		                     & -\frac{a^2k}{2h^2}(2hu_x-2(-h)^2u_{xx}/2+(-2h)u_x+(-2h)^2u_{xx}/2+O(h^3))                      \\
		                     & =u_t+\frac{1}{2}ku_{tt}+O(k^2)+au_x+O(h^2)-\frac{a^2k}{2}u_{xx}+O(h^2)                         \\
		(u_{t}=-au_{x})\quad & =O(k^2+h^2).
	\end{align*}
	For $a<0$, the proof is similar.
\end{proof}

\section{11.27 Check the stable condition for Beam-Warming and reproduce the plots.}
\begin{proof}[Proof]
	% shang dang le!
	% We only prove the case of $a>0$ . First, since a 3-step method is applied here, for convenience we suppose that, the first three $U^1,U^2,U^3$ 
	% have already been produced by other methods prior to applying our Beam-Warming method, which means the last three elements are mapped to themselves. 
	% Then observe the matrix form $U^{n+1}=AU^{n}$ of (11.28).
	% In fact our $A$ is an upper-triangular matrix and the main diagonal elements are $x=1-\frac{3\mu}{2}+\frac{\mu^2}{2}$ (except last three elements as 1 ).
	% To ensure the stability of our method, $|x|$ should be controlled within 1. Solving $|1-\frac{3\mu}{2}+\frac{\mu^2}{2}|\leq 1 (\mu\geq0)$ leads to $0\leq\mu\leq2$.

	We only prove the case of $a>0$. The matrix form $U^{n+1}=AU^{n}$ of (11.28) can be written as
	\begin{align*}
		A=\begin{bmatrix}x&y&z&&&\\&x&y&z&&\\&&\ddots&\ddots&\ddots&\\&&&x&y&z\\z&&&&x&y\\y&z&&&&x\end{bmatrix}_{(m+1)*(m+1)},
	\end{align*}
	with
	\begin{align*}
		x & =1-\frac{3\mu}{2}+\frac{\mu^2}{2}, \\
		y & =\mu(\mu-2),                       \\
		z & =\frac{1}{2}\mu(\mu-1).
	\end{align*}
	To ensure the stability of our method, the magnitude of all eigenvalues
	from A should be controlled within 1, which equals to $||A||_2\leq1$.
	Divide A into matrix $B=D+U$ and matrix $C=L$ with D, L, U being the diagonal,
	lower triangular, and upper triangular part of A. Then we have that B is upper-triangular with
	$||B||_2=|x|$ and C is lower triangular with $||C||_2=0$. As a result, $||A||_2=||B+C||_2\leq||B||_2+||C||_2=|x|$.
	Solving $|x|=|1-\frac{3\mu}{2}+\frac{\mu^2}{2}|\leq 1 (\mu\geq0)$ leads to $0\leq\mu\leq2$, which completes
	the proof.

	The case of $a<0$ is almost same and it will lead to $-2\leq\mu\leq0$.

	The plots are reproduced by MATLAB R2019b.
	\lstset{language=Matlab}
	\begin{lstlisting}
%This is a file named pro_draw_1127.m.
	function  pro_draw_1127(mu,A)
	string=['BeamWarming__mu=',num2str(mu),'.png'];
	figure(1);
	k=[0:0.01:2*pi];
	z=exp(k*i)-1.0;
	fp0=plot(z);
	fp0.Color = 'r';
	grid on;
	set(gca,'XTick',[-1.5:0.5:1.5]);   
	set(gca,'yTick',[-1.5:0.5:1.5]); 
	axis equal;
	hold on;
	E=eig(A);
	fp1=plot(E);
	fp1.Color='k';
	fp1.LineStyle='none';
	fp1.Marker='.';
	title(string);
	saveas(1,string);
	clf(1);
	end

%This is a file named pro_create_A_1127.m.
	function A = pro_create_A_1127(m,mu)
	%PRO_CREATE_A_1127 generarte A with given m and \mu.
	A=zeros(m+1,m+1);
	x=1-1.5*mu+0.5*mu*mu;
	y=2*mu-mu*mu;
	z=-0.5*mu+0.5*mu*mu;
	for i=[1:m-1]
	A(i,i)=x-1;
	A(i,i+1)=y;
	A(i,i+2)=z;
	end
	A(m,1)=z;
	A(m,m)=x-1;
	A(m,m+1)=y;
	A(m+1,1)=y;
	A(m+1,2)=z;
	A(m+1,m+1)=x-1;
	end

%This is a file named pro_11_27.m.
	m=127;
	test_mu=[0.8,1.6,2.0,2.4];
	for i=1:4
		pro_draw_1127(test_mu(i),pro_create_A_1127(m,test_mu(i)));
	end

\end{lstlisting}
	The output is shown below.

	\includegraphics[width=0.45\textwidth]{BeamWarming__mu=0.8.png}
	\includegraphics[width=0.45\textwidth]{BeamWarming__mu=1.6.png}

	\includegraphics[width=0.45\textwidth]{BeamWarming__mu=2.png}
	\includegraphics[width=0.45\textwidth]{BeamWarming__mu=2.4.png}

\end{proof}

\section{11.35 Reproduce all results in Example 11.34.}
\begin{proof}[Proof]
	To solve this initial-value advection problem, two methods can be adopted. The first way is the 'standard' way that requires determining
	the domain of dependence for the solution when T = 17 , $x \in [15,25]$. The other way, is to view this problem as 'periodic' when $x\in[-2,25]$, $t\in[0,17]$,
	and also works well. Without loss of generality, I choose the first 'standard' method.
	The plots are reproduced by MATLAB R2019b.
	\lstset{language=Matlab}
	\begin{lstlisting}
%This is a file named f0.m.	
function u = f0(x)
u=exp(-20.*(x-2).*(x-2))+exp(-(x-5).*(x-5));
end

%This is a file named f17.m.
function u = f17(x)
u=exp(-20.*(x-2-17).*(x-2-17))+exp(-(x-5-17).*(x-5-17));
end

%%The following are matlab files that generarte plots.
%initial-condition
clear
k=15:0.05:(25);
z=f17(k);
u_temp=z;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
hold on;
fp0=plot(k,z);
fp0.Color = 'none';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='b';
fp1.Marker='.';
string='Initial-condition.png';
title(string);
saveas(1,string);
clf(1);

%mu=0.8 lax-friedrichs
clear
diatah=17*5/4.0;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
n=(10+2*diatah)/0.05+1;
for i=1:25*17
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= 0.1*u(j+2)+0.9*u(j);
    end
    u=u_temp;
end
u_temp=u;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
hold on;
k=15:0.05:(25);
z=f17(k);
fp0=plot(k,z);
fp0.Color = 'b';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
string='LaxFriedrichs-mu=0.8.png';
title(string);
saveas(1,string);
clf(1);


%mu=0.8 lax-wendroff
clear
diatah=17*5/4.0;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
n=(10+2*diatah)/0.05+1;
for i=1:25*17
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= (-0.4+0.32)*u(j+2)+(0.4+0.32)*u(j)+(1-0.64)*u(j+1);
    end
    u=u_temp;
end
u_temp=u;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
hold on;
k=15:0.05:(25);
z=f17(k);
fp0=plot(k,z);
fp0.Color = 'b';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
string='LaxWendroff-mu=0.8.png';
title(string);
saveas(1,string);
clf(1);

%mu=0.8 leap-frog
clear
diatah=17*5/4.0;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
u0=u;
n=(10+2*diatah)/0.05+1;
for i=1:1
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= 0.1*u(j+2)+0.9*u(j);
    end
    u1=u_temp;
end
for i=2:25*17
    n=n-2;
    u_temp2=zeros(1,n);    
    for j=1:n
       u_temp2(j)= u0(j+2)-0.8*(u1(j+2)-u1(j));
    end
    u0=zeros(1,n);
    u0=u1;
    u1=u_temp2;
end
u_temp=u1;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
k=15:0.05:(25);
z=f17(k);
fp0=plot(k,z);
fp0.Color = 'b';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
string='Leapfrog-mu=0.8.png';
title(string);
saveas(1,string);
clf(1);

%mu=0.8 Upwind
clear
diatah=17*5/4.0;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
n=(10+2*diatah)/0.05+1;
for i=1:25*17
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= 0.2*u(j+1)+0.8*u(j);
    end
    u=u_temp;
end
u_temp=u;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
hold on;
k=15:0.05:(25);
z=f17(k);
fp0=plot(k,z);
fp0.Color = 'b';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
string='Upwind-mu=0.8.png';
title(string);
saveas(1,string);
clf(1);

%mu=0.8 Beam-warming
clear
diatah=17*5/4.0;
a=(15-2*diatah):0.05:(25);
u=f0(a);
n=(10+2*diatah)/0.05+1;
for i=1:25*17
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= (-0.4+0.32)*u(j)+(4*0.4-2*0.32)*u(j+1)+(1-0.4*3+0.32)*u(j+2);
    end
    u=u_temp;
end
u_temp=u;
figure(1);
axis([15 25 -0.4 1.2]);
set(gca,'XTick',15:1:25);   
set(gca,'yTick',-0.4:0.2:1.2); 
hold on;
k=15:0.05:(25);
z=f17(k);
fp0=plot(k,z);
fp0.Color = 'b';
grid on;
hold on;
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
string='BeamWarming-mu=0.8.png';
title(string);
saveas(1,string);
clf(1);

%mu=1 lax-wendroff
clear
diatah=17;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
n=(10+2*diatah)/0.05+1;
for i=1:20*17
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= (-0.5+0.5)*u(j+2)+(0.5+0.5)*u(j)+(1-1)*u(j+1);
    end
    u=u_temp;
end
u_temp=u;
figure(1);
set(gca,'XTick',[15:1:25]);   
set(gca,'yTick',[-0.4:0.2:1.2]); 
k=15:0.05:(25);
z=f17(k);
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
grid on;
hold on;
fp0=plot(k,z);
fp0.Color = 'b';
string='LaxWendroff-mu=1.0.png';
title(string);
saveas(1,string);
clf(1);

%mu=1 leap-frog
clear
diatah=17;
a=(15-diatah):0.05:(25+diatah);
u=f0(a);
u0=u;
n=(10+2*diatah)/0.05+1;
for i=1:1
    n=n-2;
    u_temp=zeros(1,n);
    for j=1:n
       u_temp(j)= 0.1*u(j+2)+0.9*u(j);
    end
    u1=u_temp;
end
for i=2:20*17
    n=n-2;
    u_temp2=zeros(1,n);    
    for j=1:n
       u_temp2(j)= u0(j+2)-1.0*(u1(j+2)-u1(j));
    end
    u0=zeros(1,n);
    u0=u1;
    u1=u_temp2;
end
u_temp=u1;
figure(1);
set(gca,'XTick',[15:1:25]);   
set(gca,'yTick',[-0.4:0.2:1.2]); 
k=15:0.05:(25);
z=f17(k);
fp1=plot(k,u_temp);
fp1.LineStyle='-';
fp1.Color='r';
fp1.Marker='.';
grid on;
hold on;
fp0=plot(k,z);
fp0.Color = 'b';
string='Leapfrog-mu=1.0.png';
title(string);
saveas(1,string);
clf(1);
		
\end{lstlisting}
	The output plots are shown below.

	\includegraphics[width=0.4\textwidth]{Initial-condition.png}
	\includegraphics[width=0.4\textwidth]{Leapfrog-mu=0.8.png}

	\includegraphics[width=0.4\textwidth]{LaxFriedrichs-mu=0.8.png}
	\includegraphics[width=0.4\textwidth]{LaxWendroff-mu=0.8.png}

	\includegraphics[width=0.4\textwidth]{Upwind-mu=0.8.png}
	\includegraphics[width=0.4\textwidth]{BeamWarming-mu=0.8.png}

	\includegraphics[width=0.4\textwidth]{LaxWendroff-mu=1.0.png}
	\includegraphics[width=0.4\textwidth]{Leapfrog-mu=1.0.png}

\end{proof}

\section{11.38 Derive modified equation for Lax-Wendroff.}
\begin{proof}[Proof]
	From Lemma 11.18, neglecting $O(h^4+k^4)$ term, Lax-Wendroff method can be written as
	\begin{align*}
		0        & =u_t+\frac{k}{2}u_{tt}+\frac{k^2}{6}u_{ttt}+\frac{k^3}{24}u_{tttt}+au_x+a\frac{h^2}{6}u_{xxx}-\frac{ka^2}{2}(u_{xx}+\frac{h^2}{12}u_{xxxx}), \\
		u_t+au_x & =-(\frac{k}{2}u_{tt}+\frac{k^2}{6}u_{ttt}+\frac{k^3}{24}u_{tttt})-a\frac{h^2}{6}u_{xxx}+\frac{ka^2}{2}(u_{xx}+\frac{h^2}{12}u_{xxxx}).       \\
	\end{align*}
	Since
	\begin{align*}
		u_{tt} & =-au_{xt}-(\frac{k}{2}u_{ttt}+\frac{k^2}{6}u_{tttt}+\frac{k^3}{24}u_{ttttt})-a\frac{h^2}{6}u_{xxxt}+\frac{ka^2}{2}(u_{xxt}+\frac{h^2}{12}u_{xxxxt}), \\
		u_{tx} & =-au_{xx}-(\frac{k}{2}u_{ttx}+\frac{k^2}{6}u_{tttx}+\frac{k^3}{24}u_{ttttx})-a\frac{h^2}{6}u_{xxxx}+\frac{ka^2}{2}(u_{xxx}+\frac{h^2}{12}u_{xxxxx}), \\
	\end{align*}
	we have
	\begin{align*}
		u_{tt}  & =a^2u_{xx}-a^3\frac{k}{2}u_{xxx}+\frac{ka^2}{2}u_{xxx}+\frac{ka^3}{2}u_{xxx}-\frac{ka^2}{2}u_{xxx}+O(h^2+k^2), \\
		        & =a^2u_{xx}+O(h^2+k^2),                                                                                         \\
		u_{ttt} & =-a^3u_{xxx}+O(k+h).
	\end{align*}
	Combining them together, we have
	\begin{align*}
		u_t+au_x & =-(\frac{k}{2}a^2u_{xx}+\frac{k^2}{6}(-a^3u_{xxx}))-a\frac{h^2}{6}u_{xxx}+\frac{ka^2}{2}u_{xx}+O(h^3+k^3) \\
		         & =-\frac{ah^2}{6}(1-\mu^2)u_{xxx}+O(h^3+h^3).
	\end{align*}
	Neglect $O(h^3+k^3)$ term and then we complete the proof.
\end{proof}

\section{11.40 Derive modified equation for Leapfrog and for Lax-Wendroff with one more term of higher-order derivative.}
\begin{proof}[Proof]
	Leap-frog method can be written as
	\begin{align*}
		0        & =u_t+\frac{k^2}{6}u_{ttt}+a(u_x+\frac{h^2}{6}u_{xxx})+O(h^4+k^4), \\
		u_t+au_x & =-\frac{k^2}{6}u_{ttt}-a\frac{h^2}{6}u_{xxx}+O(h^4+k^4).
	\end{align*}
	Since
	\begin{align*}
		u_{ttt} & =-au_{xtt}-\frac{k^2}{6}u_{ttttt}-a\frac{h^2}{6}u_{xxxtt}+O(h^4+k^4), \\
		u_{ttx} & =-au_{xxt}-\frac{k^2}{6}u_{ttttx}-a\frac{h^2}{6}u_{xxxxt}+O(h^4+k^4), \\
		u_{ttt} & =-au_{xxx}-\frac{k^2}{6}u_{tttxx}-a\frac{h^2}{6}u_{xxxxx}+O(h^4+k^4),
	\end{align*}
	we have
	\begin{align*}
		u_{ttt}=-a^3u_{xxx}+O(k^2+h^2).
	\end{align*}
	Combining them together, we have
	\begin{align*}
		u_t+au_x & =a^3\frac{k^2}{6}u_{xxx}-a\frac{h^2}{6}u_{ttt}+O(h^4+k^4) \\
		         & =-\frac{ah^2}{6}(1-\mu^2)u_{xxx}+O(h^4+h^4).
	\end{align*}
	Neglect $O(h^4+k^4)$ term and then we complete the proof.

	Furthermore, although it's similar to Lax-Wendroff 's modified equation, if we retain one more
	term of higher-order derivative, they will differ in the high-order term of remainder.

	For Lax-Wendroff,
	\[u_t+au_x=-\frac{ah^2}{6}(1-\mu^2)u_{xxx}+\varepsilon_fu_{xxxx},\] with $\varepsilon_f=O(k^3+h^3)$.

	For Leap-frog, however, from the symmetric form in both x and t, we see that all even-order derivatives drop out.
	Hence,
	\[u_t+au_x=-\frac{ah^2}{6}(1-\mu^2)u_{xxx}+\varepsilon_wu_{xxxxx}, \] with  $\varepsilon_w=O(k^4+h^4)$.
\end{proof}

\section{11.41 Derive modified equation for Beam-Warming and connect this to Example 11.34.}
\begin{proof}[Proof]
	From Exercise 11.26, we have
	\begin{align*}
		0        & =u_t+\frac{k}{2}u_{tt}+\frac{k^2}{6}u_{ttt}+a(u_x-\frac{h^2}{3}u_{xxx})-\frac{a^2k}{2}(u_{xx}-hu_{xxx})+O(h^3+k^3), \\
		u_t+au_x & =-\frac{k}{2}(u_{tt}+\frac{k}{3}u_{ttt})+\frac{a}{3}h^2u_{xxx}+\frac{a^2k}{2}(u_{xx}-hu_{xxx}).
	\end{align*}
	Since
	\begin{align*}
		u_{tt} & =-au_{xt}-(\frac{k}{2}u_{ttt}+\frac{k^2}{6}u_{tttt})+a\frac{h^2}{3}u_{xxxt}+\frac{ka^2}{2}(u_{xxt}-hu_{xxxxt}), \\
		u_{tx} & =-au_{xx}-(\frac{k}{2}u_{ttx}+\frac{k^2}{6}u_{tttx})+a\frac{h^2}{3}u_{xxxx}+\frac{ka^2}{2}(u_{xxx}-hu_{xxxxx}), \\
	\end{align*}
	we have
	\begin{align*}
		u_{tt}  & =a^2u_{xx}+O(h^2+k^2), \\
		u_{ttt} & =-a^3u_{xxx}+O(k+h).
	\end{align*}
	Combining them together, we have
	\begin{align*}
		u_t+au_x & =-\frac{k}{2}(a^2u_{xx}-\frac{k}{3}a^3u_{xxx})+\frac{a}{3}h^2u_{xxx}+\frac{a^2k}{2}(u_{xx}-hu_{xxx})+O(h^3+k^3). \\
		         & =\frac{ah^2}{6}(\mu^2-3\mu+2)u_{xxx}+O(h^3+h^3).
	\end{align*}
	Neglect $O(h^3+k^3)$ term and then we complete the proof.
	As a result, by Example D.27 we have
	\begin{align*}
		C_p(\xi) & =a+\frac{ah^2}{6}(\mu-1)(\mu-2)\xi^2, \\
		C_g(\xi) & =a+\frac{ah^2}{2}(\mu-1)(\mu-2)\xi^2.
	\end{align*}
	When $\mu=0.8(a>0)$, both $C_p$ and $C_g$ have a magnitude larger than a, thus
	numerical oscillations move ahead of the true wave crest; this answers Question(e).
\end{proof}

\section{11.42 Discuss $\mu=1$ for Lax-Wendroff and Leapfrog and connect this to Example 11.34.}
\begin{proof}[Proof]
	For both two methods, when $\mu=0.8$, $u_{xxx}$ term will contribute to dispersive behavior. This 
	leads to an oscillating solution, along with a shift in the location of the solution peak.
	When $\mu=1$, however, the dispersive coefficient of both methods will vanish, thus dispersive behavior will be greatly
	weakened. This answers Question(f) and explains why $\mu=1$ is much better than $\mu=0.8$ for both methods.
\end{proof}
\section{11.43 Derive ampliﬁcation factor for Upwind and show its stable region.}
\begin{proof}[Proof]
	\begin{align*}
		U^{n+1}_j & =U^n_j(1-\mu+\mu e^{-i\xi h}), \\
		g(\xi)    & =(1-\mu)+\mu e^{-i\xi h}.
	\end{align*}
	Hence, the stable region is
	\begin{align*}
		|1-2\mu| & \leq 1 ,         \\
		|\mu|    & \leq 1.(\mu > 0)
	\end{align*}
\end{proof}
\section{11.44 Derive ampliﬁcation factor for Lax-Friedrichs and show its stable region.}
\begin{proof}[Proof]
	\begin{align*}
		U^{n+1}_j & =U^n_j(\frac{1-\mu}{2}e^{i\xi h}+\frac{1+\mu}{2}e^{-i\xi h}),             \\
		g(\xi)    & =\frac{e^{i\xi h}+e^{-i\xi h}}{2}-\frac{\mu(e^{i \xi h}-e^{-i \xi h})}{2} \\
		          & =cos(\xi h)-\mu i sin (\xi h).
	\end{align*}
	Hence, the stable region is
	\begin{align*}
		|\mu|\leq 1.
	\end{align*}
\end{proof}
\section{11.45 Derive ampliﬁcation factor for Lax-Wendroff and show its stable region.}
\begin{proof}[Proof]
	\begin{align*}
		U^{n+1}_j & =U^n_j(1-\mu^2-\frac{\mu-\mu^2}{2}e^{i\xi h}+\frac{\mu^2+\mu}{2}e^{-i\xi h}),          \\
		g(\xi)    & =1-\mu^2-\frac{\mu}{2}(e^{i\xi h}-e^{-i\xi h})+\frac{\mu^2}{2}(e^{i\xi h}+e^{-i\xi h}) \\
		          & =1-\mu^2+\mu^2cos(\xi h)-\mu sin(\xi h)i                                               \\
		          & =1-2\mu^2sin^2\frac{\xi h}{2} - i\mu sin(\xi h).
	\end{align*}
	Solve the inequality, and the stable region is
	\begin{align*}
		|\mu|\leq 1.
	\end{align*}
\end{proof}
\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
